Spin correlations and doublon production rate for fermionic atoms in modulated 
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We compute the integrated doublon production rate in response to a lattice modulation for two- 
component fermions in an optical lattice. We derive a general formula for the integrated intensity, 
valid in the presence of inhomogeneous potentials such as the trap, which gives the integrated 
intensity in terms of equal time correlation functions only. Such a formula is thus well suited for 
direct numerical calculations. We show that, in the limit of large repulsion for commensurate fillings, 
or for temperature ranges for which the hopping is incoherent, the integrated doublon spectrum is 
directly related to the nearest-neighbor spin-spin correlation function. We compute its temperature 
dependence in this regime using finite temperature quantum Monte Carlo calculation. 
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Ultra-cold atom gases confined in an optical lattice 
provide a controlled realization of the Hubbard model pj, 
|2| which plays a central role in the study of strongly cor- 
related electron systems such as Mott transition, high-Tc 
superconductivity, and quantum magnetism. Despite in- 
tensive research the phase diagram of this model is still 
largely debated. Cold atoms, with the control of the in- 
teractions using a Feshbach resonance technique and the 
control of kinetic energy via lattice depth, allow us to 
probe such physics in detail. 

In order to do so, the development of experimental 
probes to capture many-body quantum states is also 
an important issue. In connection with the Hubbard 
physics, measuring the antiferromagnetic (AFM) corre- 
lations is of prime importance. However this is not an 
easy task. The AFM order can potentially be extracted 
from time-of-flight measurements via shot-noise measure- 
ment [H-IH. However, the measurement is far from trivial 
and also depends crucially on the direction of the mag- 
netic order. Local addressing [1, 0] is also a potential 
route but systems are for the moment quite small, and 
in addition such a technique is complicated to extend to 
three-dimensional systems. It is thus highly desirable to 
study other probes that can potentially give direct access 
to the magnetic correlations 

One probe which has proved to be very efficient and 
relatively simple to implement is the amplitude modula- 
tion of the optical lattice fis'l. This probe, first used for 
bosonic systems with a measurement of the absorbed en- 
ergy, gives direct access to the kinetic energy correlation 
functions [l6l - [T8| . For fermions the energy absorption 
rate (EAR) cannot be implemented accurately enough 
and a variant of this probe measuring the doublon pro- 
duction rate (DPR) has been proposed [l^. This DPR 
measurement in the linear response regime could be suc- 
cessfully implemented ^2^, and thus stimulates further 
studies of theoretical calculation of DPR spectra (2l| - [26| 
and doublon dynamics p7l - [29| . 

In addition, it was shown [Toj by a direct compari- 
son of the two quantities, that the integrated intensity of 



the modulation does directly give access to the nearest- 
neighbor AFM correlations. Intuitively, for fermionic 
atoms, spin configurations of neighboring atoms are rel- 
evant to the spectra because of the Pauli exclusive prin- 
ciple: while the hopping is allowed for neighboring atom 
spins pointing in anti-parallel, it is blocked for ferromag- 
netically aligned spin configuration. This point, which 
allows one to use the modulation as a simple probe of 
magnetism, was explored further both theoretically and 
experimentally in the linear response regime at high tem- 
peratures [30| . 



In this paper, we further analyze the integrated in- 
tensity of the DPR. We consider potentially inhomoge- 
neous systems, e.g. due to the presence of the trap. We 
show that the frequency integrated DPR can be fully ex- 
pressed, within linear response, in terms of static correla- 
tions. Such correlations, contrarily to the original DPR 
response at fixed frequency, are well within the reach, 
without any need for analytical continuation, of pow- 
erful numerical methods such as quantum Monte Carlo 
(QMC) simulation, thus potentially allowing a very pre- 
cise comparison of the integrated DPR and theoretical 
calculations. For the case when the filling is commen- 
surate and the interactions are large compared to the 
kinetic energy, or when the hopping is incoherent due 
to the temperature, that the above formulas reduce, in 
agreement with the initial study of Ref. [l9|, to a mea- 
sure of the nearest-neighbor AFM correlation functions. 
We obtain the temperature dependence of this quantity 
via a QMC simulation of the Heisenberg quantum spin 
model in a three-dimensional cubic lattice. 



We consider two-component fermionic atoms strongly 
confined in an optical lattice. In the case of deep opti- 
cal lattices, the physics is well described by the spin- 1/2 
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fermionic Hubbard model, Hq = Hk + -ffu 



_ffp, with 



(1) 
(2) 
(3) 



where Cja and are, respectively, an annihilation and 
number operator of a fermionic atom at a jth site. The 
on-site potential Vj (for example, corresponding to trap 
potential) generally breaks translational symmetry. 

The parameters J and U are given as a function 
of an amplitude of a sinusoidal lattice potential, i.e., 
K)p (^") ~ Vb [cos^ {txx /a) + cos^ (Try /a) + cos^ (ttz /a)] , 
where the lattice constant is a [IJ. The dynamical am- 
plitude modulation defined on continuum space is rep- 
resented in an effective lattice model as follows: J ^- 
J[l-|-(5Jcos(wi)] and U — >■ U[l + 5U cos{ujt)\ correspond- 
ing to Vb — )■ Vb -I- 5V cosiujt) for the lattice depth [3l|. 
Thus the perturbation Hamiltonian of the lattice modula- 
tion is given as V(t) = {6 J) cos{ujt)HK + {SU) cos{ujt)Hij, 
but here we use the alternative representation to simplify 
the problem [l^ as V{t) = (SU) cos(iui)i7o + cos{ujt)S 
with 



S ^ idF)HK ~ {SU)Hp, 



(4) 



where SF = SJ — SU is a dimensionless perturbation 
parameter. Note that in the homogeneous case (vj = 
v) the perturbation operator becomes the kinetic energy 
Hk. 

In the linear response regime, the DPR per site can 
be defined as an increment of the doublon number time- 
averaged over a single period 27r/a;, 
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f7 27r/cj 



*+(2^/^) ^ 

dt' ^Nu{t'), (5) 



where A^D(i) is the number of created doublons, and 
ri is a total site number of the system. Note that as 
seen below, the time dependency in the right-hand side 
should cancel due to the single-period time average in 
the linear response region. The doublon number can 
also be written by the Hubbard interaction A^d(^) — 
{H\j) /U, which is averaged by the density matrix p{t) — 
e-HW/ksT ^rj.^^^-H{t)/kBT^ for time-dependent Hamilto- 
nian H{t) = Ho + V{t). Thus the doublon number is also 
expressed as 



E{t)-{HK+H^)-{V{t)) 
U 



(6) 



where E(t) = {H{t)) is a system energy. In the second- 
order perturbation theory with respect to the lattice 
modulation V{t), all terms in Eq. ^ apart from the first 
one are found to generate only oscillatory terms whose 



contribution to the DPR spectrum disappears due to the 
time average in Eq. ([5]) . Therefore one can rewrite Eq. ([S]) 
as follows: 
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t+(27r/Lj) 



dt' 
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-E{t'). (7) 



Interestingly, the equality means that the DPR is equiv- 
alent to the EAR in the second-order response regime, as 
was pointed out in Ref. [19]. The EAR ([7]) as a second- 
order response can be formulated by the linear response 
theory. As a result, the DPR formula can be given [f^ 
as 



ft(c.) = -^^Im[.xtM] 



(8) 



where = -i dt e"^* {[S{t), S{0)])q is the Fourier 

transform of the retarded correlation function of the op- 
erator S for the unperturbed Hamiltonian Hq. 

We now consider the general formula of DPR spectra 
([8]) integrated over the modulation frequency uj. Before 
implementing the integral, we clarify the point of the in- 
tegral range of modulation frequency from an experimen- 
tal point of view. In actual experiments a high-frequency 
cutoff is necessary. In order to discuss the needed cut- 
off, let us recall the setup of the Hubbard model. Atoms 
in an optical lattice form a Bloch band structure, and 
the band gaps are determined by depth of the optical 
lattice potential. Since the single-band Hubbard model 
Hq is introduced to demonstrate the physics in the low- 
est Bloch band, the high energy cutoff necessary to jus- 
tify the effective model should be taken to be sufficiently 
small compared with the band gap between the lowest 
and the next Bloch band. Therefore, by integrating the 
DPR spectrum over the frequency region below the band 
gap, the integration discussed below can be estimated in 
experiments. 

The integrated DPR (O reads 



duj uJXsi^) 



dt 



{[[Ho,S{t)],Sm)o 
t + iO+ 



(9) 



Here the quantity ([[i/o, S'(i)], 5(0)])o in Eq. ^ is real. 
Thus the imaginary part in the right-hand side of Eq. Q 
is caused by the denominator of the integrand. Using the 
decomposition, l/(t-t-iO+) = Vj—iTrS{t) where 7-" denotes 
the principal value, one can take the imaginary part, and 
obtain the dimensionless total spectrum weight T as 



F = 



did 



iJ/h) 



1 TT{[[Hn,S],S])o 

2 ^D^'^) - ^ • (10) 



This formula is one of the central results of this paper. 
Remarkably, Eq. ([T0|) is deduced from Eq. ([8]) without 
any approximation, and in addition the inhomogeneity 
effect of the system is fully taken into account in the for- 
mula. It relates a dynamical quantity to a calculation of 
equilibrium correlation functions. Such correlations are 
easily amenable to a computation via a large number of 
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numerical techniques such as QMC simulations, allowing 
us to make direct contact between the measured DPR 
spectrum and various quantities of the Hubbard model. 
One can thus, by comparing the experimental and theo- 
retical results, expect to determine hard to get parame- 
ters, such as the temperature (or entropy). 

Furthermore, as we will show below this general cor- 
relation, the generic formula (jlO[) is reduced in the inter- 
esting regime of temperatures fc^T -C ?7 to a direct mea- 
sure of the spin-spin correlations. Equation ((TU]) contains 
two-particle correlation functions such as (c|^c^j^ 0^^04)0. 
Thus the practical calculation is not trivial in general, 
and the physical interpretation is not also clear. In what 
follows, we separately consider the two limited but inter- 
esting regions: (i) high-temperature /cbT ^ J and (ii) 
strongly correlated regime J <^ U aX half-filling. Note 
that the two-particle correlation functions can then be 
written as 

+ 5i,i5j^k{c\^c]^Cj^c^{)Q. (11) 

In the first case (i), hopping is incoherent regardless 
of the filling. Thus we can simply apply Eq. (|lll) to the 
generic formula ()10|) . Consequently F simplifies as 



2'K{5Ff 
t:{5F){5J) 



where n'; 
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V, ~ Vj 



U 



(K)o - (<)o) 



(12) 



and 2; is a coordination number. 



The total doublon number in equilibrium state, Nd ~ 
^j{nj^nj^)Q, has also been introduced. Note that the 
inhomogeneity of a system is also fully taken into account 
in formula (jl2p . In the homogeneous case, since Vi — Vj ~ 
0, F turns out to be given only by the first line. 

In the second strongly interacting regime (ii), which 
allows us to reach much lower temperatures, one has to 
restrict the study to the half-filling homogeneous case to 
avoid regions in which coherent hopping could still exist. 
In particular, in the presence of a trap, this will occur 
in the shell of the compressible region. Hopefully such 
regions would give small contributions compared to the 
bulk of the response. In addition, the response of these 
regions is mostly concentrated at low energy rather than 
for energy or order U, and thus most of their contribu- 
tion can be filtered away in the frequency integration. 
If we restrict to the commensurate case, the simplifica- 
tion (ITTT) is still applicable for ksT ^ U. Then the on- 
site potential term Hp essentially disappears (3^ . and 
the integrated spectrum F turns out to be Eq. ([TU]) in 
which S is replaced by {SF)Hk- The resultant F after 
applying Eq. (in]) is given by the first line of Eq. (IT^ . 
However, the density fluctuation would then be strongly 
suppressed, and one can approximately take Nd w and 



(njnpo ~ 1, whose temperature dependency would be 
negligibly small. Therefore, the integrated DPR spec- 
trum is given as 



2TTz{8Ff 



(13) 



where z is a coordination number. What is important 
is that only the nearest-neighbor spin-correlation func- 
tion {Si ■ Sj) is a dominant contribution to the tempera- 
ture dependence of F. Namely, F can be identical to the 
nearest-neighbor spin correlation. 

The nearest-neighbor spin correlation appearing in 
Eq. (fT5)) is exactly the same as the energy of the quan- 
tum spin Heisenberg model. As is well known, in the 
strongly interacting region fc^T, J ^ U, at half-filling, 
the Heisenberg model is deduced from the Hubbard 
model as a consequence of the second-order perturbation 
theory in terms oi J/U . Thus, in such a region, Eq. (|13p 
allows us to probe the system energy. 

In order to see the temperature dependence of F, we 
numerically estimate the temperature dependence of the 
integrated DPR spectrum F. However, it is a highly non- 
trivial problem because it is necessary to compute the 
dynamical quantity (|S]) directly in a wide regime of tem- 
perature. Thus, in this paper the parameter region is re- 
stricted to the regime for which the approximation (jl3p 
is valid, i.e., fc^T, J <C J7. In addition, we calculate 
the static quantity (nearest-neighbor spin correlation) 
{Si ■ Sj) instead of integrating the DPR spectrum over 
the modulation frequency. In the low-temperature Mott 
insulating case, we can approximately choose the AFM 
Heisenberg model as an effective model, which allows one 
to compute the nearest-neighbor spin correlation through 
the energy of the Heisenberg model. 

We implement the QMC calculation for the AFM 
Heisenberg model of a cubic lattice where the system 
size is 17 = 14'^ [S^]. The result is shown in Fig. [H where 
the calculated energy is translated into F by use of the 
relation (|13p . As a reference of the critical temperature 
of the Neel transition [s^, [1^ , the specific heat as a func- 
tion of temperature is also shown. The derivative of F 
as temperature rises is found to be maximum near the 
critical point, which seems an inflection point. On the 
other hand, the nearest-neighbor spin correlation and the 
energy for temperatures higher than the critical point ex- 
hibit long tail decay. 

To summarize, we have explicitly discussed the DPR 
spectra in the presence of an inhomogeneous potential, 
and the generalized form of the DPR spectrum as a func- 
tion of lattice modulation frequency has been obtained. 
In addition we have implemented the integration of the 
DPR spectrum over modulation frequency. The problem 
of dynamics can be reduced to the calculation of static 
quantity as shown in Eq. (jlOp . Such a formula can be the 
basis of a very accurate numerical comparison with the 
experiment. Furthermore, focusing on the two different 
regimes, (i) fc^T ^ J and (ii) J ^ ?7, a generic form of 
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the integrated DPR spectrum has been deduced. In par- 
ticular, for the incompressible Mott regime, it is found to 
be associated with the nearest-neighbor spin correlation 
as in Eq. ([T^. and can thus be used as a probe of the 
magnetic properties. In the larger J regime, higher order 
hopping processes would be more relevant, and the inte- 
grated DPR spectrum is expected to probe longer-range 
spin correlations. One can thus use the shaking of the 
lattice and DPR to directly probe for magnetic order. 



FIG. 1. (Color online) The temperature dependence of (top) 
the integrated DPR spectrum per site F scaled by 2-kz{SF)^, 
(middle) the energy per site scaled by a superexchange spin 
coupling Jox, which is identical to the nearest-neighbor spin 
correlation per site multiplied by the coordination number 
z = 6 in a cubic lattice, and (bottom) the specific heat. 
The QMC calculation is implemented for the quantities in 
the Heisenberg model for 14^ system size. While error bars 
are also shown, they are negligibly small. The Neel transition 
occurs near ksT ~ Jcx, and then the integrated spectrum and 
the energy per site exhibit an inflection-point-like behavior. 
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